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ABSTRACT 

We present a high resolution study of the impact of realistic satellite galaxies, extracted 
from cosmological simulations of Milky Way haloes including 6 Aquarius suites and 
Via Lactea II, on the dynamics of the galactic disc. The initial conditions for the 
multi-component Milky Way galaxy were generated using the GallC code, to ensure a 
system in dynamical equilibrium state prior to addition of satellites. Candidate sub¬ 
haloes that came closer than 25 kpc to the centre of the host DM haloes with initial 
mass enclosed within the tidal radius, M t ;d > 1O 8 M 0 =0.003 Mdi sc , were identified, 
inserted into our high resolution N-body simulations and evolved for 2 Gyr. We quan¬ 
tified the vertical heating due to such impacts by measuring the disc thickness and 
squared vertical velocity dispersion er 2 across the disc. According to our analysis the 
strength of heating is strongly dependent on the high mass end of the subhalo distri¬ 
bution from cosmological simulations. The mean increase of the vertical dispersion is 
~ 20 km 2 s -2 Gyr -1 for R > 4 kpc with a flat radial profile while, excluding Aq-F2 
results, the mean heating is < 12 km 2 s -2 Gyr -1 , corresponding to 28% and 17% of 
the observed vertical heating rate in the solar neighbourhood. Taking into account the 
statistical dispersion around the mean we miss the observed heating rate by more than 
3er. We observed a general flaring of the disc height in the case of all 7 simulations in 
the outer disc. 
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1 INTRODUCTION 


According to the A cold dark matter (ACDM) paradigm, 
structure formation follows a hierarchical scenario; dark 
matter halo growth occurs via the accretion of smaller 


size systems (e.g. White & Rees 1978 Moore et al. 19991. 
High-resolution cosmological simulations have demonstrated 
that DM substructures survive within environments sim¬ 
ilar to the Milky Way ( Moore et ak| 19991. Galaxy merg¬ 
ers are believed to be the main drivers for significant kine- 
matical perturbations of galaxies ( Toomre fc Toomre|1972 1; 
these can be classified into three categories. With regards 
to theoretical (|Hernquist|1992[) and observational evidences, 


(e.g. Woods, Geller fc Barton|[2006 Woods fc Geller|[2007 1 
“major” mergers (mass ratios > 1:3 of total galaxy mass 


* E-mail: reza@ari.uni-heidelberg.de (RM); just@ari.uni- 

heidelberg.de (AJ) 

f Fellow of the International Max Planck Research School for 
Astronomy and Cosmic Physics at the University of Heidelberg 
(IMPRS-HD). 


strongly perturb the system and can even change the mor¬ 
phology. However such events are less common compared to 
“minor” mergers with masses in the range 1:3 - 1:50 of total 
mass - 1:1 of disc mass - which are likely to destroy thin 
discs (Kazantzidis et al. 2008]). The third type, which forms 
the basis of this study, is the interaction of satellite galax¬ 
ies with intermediate mass ratios which are more commonly 
represented in units of the disc mass (< 1:10 of disc mass). 
Our Milky Way galaxy has been experiencing such merg¬ 
ers; observations of stellar tidal stream remnants provide 
evidence for such encounters (e.g. Ilbata, Gilmore & Irwin 


1994 Majewski et al.||2003 Martfnez-Delgado et al.||2005 |. 
Understanding the evolution of disc galaxies has challenged 
both observers and theorists for the past two decades. There 
have been numerous attempts at bridging the observed char¬ 
acteristics of spiral galaxies and their violent history (e.g. 
|Navarro fc White|1994| |Abadi et al.|2003|). So me of the early 
studies, including that by Moore et al. (19991 indicated that 


interaction of substructures with galactic discs could induce 
heating processes resulting in thickening of the disc. Also the 
ACDM scenario predicts a continuous infall of subhaloes - 
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including massive satellites, which are expected to impact 


the dynamics of the stellar disc significantly (Stewart et al. 


20081. However early observations of our solar neighbour¬ 
hood found the presence of a 10 Gyr old “thin” disc with 


a scale height < 500 pc (e.g. Kent et al. 1991 Dehnen & 
Binney|1998 Just fc Jahreifi|[2010 1. 


The survival of a thin disc in the (ACDM) universe 


has been tested through early analytic models (e.g. Toth 


& Ostriker 1992 

Sellwood et al. 1998) together with N- 

body/hydrodynamical simulations (e.g Quinn & Goodman 

1986; Quinn, Hernquist & Fuliagar|| 1993. 

Walker, Mihos & 

Hernquist||1986| Velazquez & Whitejl999 

Sommer-Larsen, 

Gotz & Portinari 

2003). One of the major shortcomings 


of analytical approaches was the choice of nearly circular 
satellite orbits. This is not in agreement with the ACDM 
results, where eccentric orbits are favoured (jGhigna et al. 


19981. Another disadvantage was the use of a rigid po¬ 


tential for the host halo, disc and/or bulge instead of live 
components. Rigid potentials suppress the effect of inter¬ 
actions between the halo and the disc, as well as the dy¬ 
namical friction - and hence transfer of angular momen¬ 
tum - between different components as the satellites spiral 


towards the centre. Velazquez & White (1999) have shown 


that using a rigid potential instead of a live halo overesti¬ 
mates the heating of the galactic disc by a factor of 1.5- 
2. Early numerical simulations suffered from an unrealistic 
high abundance of substructure compared to the ACDM cos¬ 
mology prediction. In addition to galaxy mergers, other pro¬ 
posed mechanisms are expected to contribute towards the 


heating. For example, Carlberg (1987) andSellwood (2013) 


looked into the induced changes of disc’s phase-space dis¬ 
tribution due to transient spiral waves and concluded that 
such structures are able to increase the velocity dispersion of 
stars, although in the disc plane, an increase of the velocity 
anisotropy and hence smaller a z /<jr is observed. Scattering 
due to Giant Molecular Clouds (GMCs) is another possible 
scenario (Kokubo & Ida 1992). An overview of possible disc 


heating mechanisms are discussed in the review by Sellwood] 
(2014). 


More recent simulations (e.g. Font et al 12001 

Ardi et 

al.J2003; Kazantzidis et al. 2008 Read et al.||2C)08 

Villalo- 

bos & Helmi|2008 

Kazantzidis et al. 2009 

1 have shown that 


infall of satellite galaxies could result in either destroying 
or greatly heating the disc, forming a thick component; this 
contradicts our picture of the Milky Way thin disc compo¬ 
nent. It is important to note these analyses focussed on the 
impact of satellites at the massive end of the substr ucture 
spectrum, with masses > 5 x 10 9 Mg. For instance, Read 


et al. |2008| concluded that accretion of a satellite with 
mass comparable to the Large Magellanic Cloud (LMC) - 
10 10 Mq - is sufficient to explain the thick disc formation. 
Also |Moster et aTj ( |2010| ) examined the role of a gas com¬ 
ponent in disc instability and found, that for a gas mass 
fraction of 25 percent, the thickening of the disc could de¬ 
crease by 20 percent - thus suppressing vertical heating. In 


addition, Quillen (2009) and Bird et al. (20121 showed that 


satellites infall could also induce radial migration in MW-like 
systems with greater impact at larger radii, while |Minchev et| 
al. (20121 andVera-Ciro et al. (20141 concluded that radial 


migration has minimal impact on the thickening of the MW 
disc. The radial migration as a result of mergers reduces 
the disc flaring (Minchev et al.|20141 which provides an op¬ 


posite picture to the measured one, especially in the outer 
parts, in the absence of external perturbations as demon¬ 


strated by Minchev et al. (20121. Based on perturbation 
theory, Vesperini & Weinberg (2000) have shown that en¬ 


ergy deposit in the perturbation by satellites, scales as the 
square of the satellite mass, and that in the resonance situ¬ 
ation, minor mergers may have a measurable impact on the 
host galaxy. In a cosmologically based hydrodynamical sim¬ 


ulation, Gomez et al. (2016) have shown that even a distant 
encounter of a satellite galaxy with M = 4 x 10 1o Mq (~ 
5 per cent of the host galaxy mass) at R = 80 kpc can per¬ 
turb the outer disc at R > 12 kpc significantly and form a 
Monoceros ring-like warp. 

The observations of the local stellar kinematics have 
suggested multiple scenarios for the evolution of the thin 
disc. The most useful observational quantity in order to in¬ 
vestigate the evolution of our Milky Way disc is the so-called 
” age-velocity dispersion” relation. This relation is aimed at 
connecting stars’ kinematics with their dynamical history 
and providing constraints on the vertical structure of the 


Galactic disc. Nordstrom et al. (2004) and follow-up studies 


(e.g. Holmberg, Nordstrom & Anderson 2007 20091 have 


favoured a picture in which the disc has experienced con¬ 
tinuous heating through its lifetime. Another possible sce¬ 
nario is the saturation of vertical heating of stars older than 
4.5 Gyr, as suggested by a detailed analysis of the Geneva- 


Copenhagen survey (Scabrokc & Gilmore 20071. Holmberg et 


al. (2009) fitted power laws to the radial, tangential, vertical 
and total velocity dispersions of ~ 2600 F and G dwarfs in 
the solar neighbourhood and derived power law exponents of 
7 = 0.39, 0.40, 0.53 and 0.40 for these respective quantities. 

Within the framework of this study, we use high- 
resolution N-body simulations (Bien et al. 20131 to investi¬ 
gate the impact of accreting satellite galaxies on the dynam¬ 
ics of our Milky Way “thin” disc and answer the question: 
are we able to explain the observed heating of the Galac¬ 
tic disc in the presence of satellites infall? For improved 
realism, the initial conditions (ICs) of our satellites follow 
that of the ACDM model ( Springel et al.|[2008 l. The inter¬ 
actions of satellites and the stellar disc have the potential 
to trigger different types of both small and global scale per¬ 
turbations, such as initiating modes in the Galactic disc. 
Hence, for a deeper understanding of such distortions in the 
phase-space distribution of disc stars, we require simulations 
that are capable of resolving the dynamics down to few tens 
of parsecs scale; a resolution much higher than most cur¬ 
rent Milky Way interaction analysis simulations. We employ 
multi-component models for our Milky Way which enable us 
to use both observational and cosmological constraints for 
modelling our Galaxy (Yurin & S pringel|20 141. 

This work improves the current knowledge about the 
impact of subhalo interactions on the observed vertical 
structure of our Galactic disc. Such structure is the result 
of a complex interplay of different heating mechanisms and, 
since we do not yet know about the exact contribution from 
each of such mechanisms, quantifying one will provide a 
more accurate picture of other processes. For the first time, 
we have analysed and compared the impact of subhaloes 
from Aquarius and Via Lactea II Milky Way-like simulation 


suites (e.g. Springel et al. 2008 Diemand et al. 20081 to 
provide a statistically valid study of satellite impact on the 
vertical structure of the stellar disc of our Galaxy. In addi- 
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Table 1. the cosmological initial conditions 


Cosmological Suite 

Om 

Da 

(78 

n s 

Aquarius 

0.25 

0.75 

0.9 

1.0 

VLII 

0.238 

0.762 

0.74 

0.951 


tion, having initial conditions for the MW multi-component 
galaxy which are in much better equilibrium state, com¬ 
pared to previous work, allows us to minimize the bias ef¬ 
fect which, in principle, could overshadow the expected real 
heating measurements. The combination of 40 pc vertical 
resolution together with ICs in a robust equilibrium state 
differentiates this analysis compared to previous studies - 
allowing minor vertical heating effects to be measured and 
distinguished from other processes. 

In section[2] we discuss the cosmological simulation 
suites that we use for the purpose of our realistic initial con¬ 
ditions of DM structures together with the statistics of such 
objects. The N-body model used for the multi-component 
Milky Way and satellites, a description of our N-body code 
and the isolated systems simulation results are explained in 
section[3] The results of simulations with satellites for the 
solar neighbourhood and across the disc, together with an 
analysis of the orientation of subhalo distribution, are dis¬ 
cussed in section[4] and we finish the paper with a summary 
in section[5] 


2 COSMOLOGICAL SIMULATION SUITES 


For a realistic picture of the interaction between satellite 
galaxies and the Milky Way’s thin disc we have extracted 
the distribution of subhaloes from two suites of cosmolog¬ 
ical simulations: Aquarius and Via Lactea II (e.g. Springel 
et al.||2008 Diemand et al.||~2008 l. The major difference be 
tween the Aquarius and Via Lactea II (hereafter VLII) lies in 
their adopted cosmology, which in the case of Aquarius are 


the WMAP one-year results (Klypin et al. 1999), while for 


VLII WMAP three-year results were used. Hence the value 
of as, which represents the amplitude of the (linear) power 
spectrum on the scale of 8 h _ 1 Mpc, is lower for VLII. This 
parameter plays a crucial role in the context of influencing 
structure growth at early epoch. Also n s corresponds to the 
spectral index of scalar fluctuations (Tablc[l]) and is close 
to being scale-invariant. The value for fl m , ag and n a are 
smaller for VLII compared to Aquarius which might impact 
the substructure distribution at the very low mass end. 


2.1 Host dark matter halo properties 

The Aquarius suite consists of a set of 6 simulations; each 
associated with a different realisation of DM halo likely to 
host a Milky Way-like galaxy and which has not experienced 
a recent major merger. The haloes are chosen from a par¬ 
ent halo of homogeneous resolution. These simulations are 
performed in cubic box with side length of 137 Mpc while 
VLII is performed in a smaller box of length 40 Mpc. The 
second highest resolution set of Aquarius simulations “Aq- 
A2- • • Aq-F2” are employed for this analysis. In our study 
we use the 2 = 0 snapshots from all the above simulations, 
which resemble the present day distribution of subhaloes. 


Each snapshot contains information for the subhaloes - in¬ 
cluding the position, velocity, maximum circular velocity 
Vmax, position of this velocity r vmax and tidal mass M t id 
for every substructure residing within the simulation box. 

The parent haloes from these seven cosmological simu¬ 
lations do not possess exactly similar profile characteristics 
such as M 200 , the mass enclosed within the sphere with ra¬ 
dius r 2 oo, where r 2 oo represents the radius at which the mean 
density of the DM halo is 200 times the critical density of 
the universe (p cr it). Thus, we scale the properties of sub¬ 
haloes from all seven simulations in terms of their position, 
velocity and mass using similar scaling relations mentioned 


by Kannan et al. (2012). All the simulations were scaled with 
respect to the M 200 of the “Aq-D 2 ” simulation - 1.77 x 10 12 
M©. The quantity g represents the scaling factor in 


M = M orig /p, 


(1) 


uorig,x,y,z 

0 V3 


and 


( 2 ) 


y : 


Vorig _ Zorig , . 

gl/i y “ gl/3 ~ “ gl/3 ' W 

The subscript “orig” corresponds to the original non-scaled 
quantities. Such scaling is important for a statistically mean¬ 
ingful comparison of different simulations. After having done 
the rescaling, all the host haloes have the same mass and size 
but different concentrations, hence different r vm ax and Knax- 
According to early N-body simulations of dark matter 
haloes following the ACDM structure formation model, the 
density profile of the DM haloes could be well described via 


the known Navarro-Frank-White (NFW) profile (Navarro et 
~ Ifl997 1 

P= 


PnfwM = 


(r/r s )(l + ryV s ) 2 ' 


( 4 ) 


where p s and r s are the scale density and scale radius of the 
subhalo. The p s can be computed given the concentration 
c of the halo which is simply the ratio of r 2 oo/r s and the 
overdensity 5 C with respect to p cr it using 


5 C = 


200 


_ _ _ (s') 

pcrit 3 ln(l + c) — c/(l + c) 

Tablc[2] shows the normalised properties of the host haloes 
from VLII and six Aquarius simulations. The particle resolu¬ 
tion of the simulation suites m p , or i g which is ~ 10 4 Mq moti¬ 
vated us to perform a mass cut and consider subhaloes with 
Mtid ^ 10 6 Mq since substructures with lower mass consist 
of only a few hundred particles. This introduces some uncer¬ 
tainties in determining the characteristics of these subhaloes 
at the outer regions. The analysis performed on the sub¬ 
haloes in this paper include this lower mass cut. In addition 
to M 200 and r 2 oo, the maximum circular velocity V ma x and 
concentration c are listed. For each parent halo we also know 
r 5 o, the radius at which the mean density is 50 times p cr it; 
V 50 , the circular velocity at this radii and f su h = / s c ub"( r so), 
the percentage fraction of the cumulative mass residing in 
substructures relative to the host halo enclosed mass at r 50 
and seen as a measure of the host halo dumpiness. It appears 
that VLII has less mass residing in bound substructures - 
/sub ~ 5 percent - compared to the Aquarius simulations. 

Fig. a illustrates / a c p ™(r), which is the substructure 
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Table 2. Normalised main halo properties; m p is the original particle mass and g is the scaling factor. M200 is the normalised mass 
enclosed within a sphere of radius r2oo, where the mean density of the halo is 200 times the critical density of the universe, c is the 
concentration while V ma . x is the maximum circular velocity of the halo, rso is the radius of the main halo at which the halo encloses the 
mass with mean density 50 times the background density and V50 is the circular velocity at this radii. / su b represents the cumulative 
mass fraction which resides in substructures relative to the total enclosed mass within rso- 


Simulation 

mp,orig 

[M 0 ] 

9 

M200 
[10 12 M 0 ] 

^200 

[kpc] 

c 

V" max 

[kms -1 ] 

^50 

[kpc] 

V50 

[kms -1 ] 

/sub 

[%] 

VLII 

4.10 x 10 3 

1.188 

1.774 

242.8 

13.29 

212.9 

425.7 

151.8 

4.77 

Aq-A2 

1.37 x 10 4 

0.963 

1.774 

242.8 

16.19 

205.9 

428.2 

156.3 

12.46 

Aq-B2 

6.45 x 10 3 

2.166 

1.774 

242.8 

9.72 

204.0 

418.0 

152.6 

9.61 

Aq-C2 

1.40 x 10 4 

1.000 

1.774 

242.8 

15.21 

222.4 

417.1 

152.3 

6.67 

Aq-D2 

1.40 x 10 3 

- 

1.774 

242.8 

9.37 

203.2 

425.7 

158.1 

13.14 

Aq-E2 

9.59 x 10 3 

1.497 

1.774 

242.8 

8.26 

204.8 

421.3 

153.8 

9.85 

Aq-F2 

6.78 x 10 3 

1.564 

1.774 

242.8 

9.79 

196.3 

424.7 

155.0 

12.80 




Figure 1. Cumulative subhalo mass fraction, substructure mass 
relative to the halo’s enclosed mass at r, as function of distance 
from the host halo centre. 


mass relative to the halo enclosed mass at r, as function 
of distance from the parent halo’s centre for all seven simu¬ 
lations. Within the inner 20kpc Aq-B2 has the highest sub¬ 
halo mass fraction, while Aq-C2 possesses the lowest value 
across nearly all radii. Fig.[2]illustrates an inverse cumulative 
histogram for the number of subhaloes as function of their 
normalised M t id- All the simulations have a similar slope 
while the scatter in the high-mass regime is the result of low 
number statistics. For VLII the slope decreases at a subhalo 
mass < 2 x 10 6 M@; below this mass-range the difference is 
sensitive to the subhalo finding algorithm and also the mass 
resolution. Accordingly, the substructure distribution is in¬ 
fluenced by lower mass subhaloes that could dominate over 
higher mass ones by few orders of magnitude in abundance. 


2.2 Identifying the most effective subhalo 
candidates 

As mentioned above, we are interested in subhaloes which 
are likely to interact with the galactic disc; hence we need to 
identify such subhaloes. We apply a simple selection crite- 


Figure 2. Inverse cumulative histogram for the number of sub¬ 
haloes as function of subhalo mass within rso °f the host haloes. 
The crossed population corresponds to the subhaloes which come 
closer than 25 kpc to the centre of the host halo within 2 Gyr 
while all represents the total sample. The vertical line marks the 
mass cut at 10® Mg. 


rion by estimating the pericentre distances of all subhaloes. 
We select these by determining the orbit of subhaloes in 
the presence of the background potential of the host DM 
halo only. Using the initial 6-dimensional phase-space in¬ 
formation of the subhaloes we integrate their orbits using 
the standard Runge-Kutta 4th order orbit integrator, which 
determines the position and velocity of subhaloes at each 
integration step (Press et al. 19961, for 2 Gyr. Any subhalo 
that comes closer than 25 kpc to the parent halo’s centre 
at any point during its orbit is marked as a crossed candi¬ 
date. This calculation does not take into account dynamical 
friction and the gravitational force of the disc and bulge. 

Fig. 0 shows Mtid vs. r per i for subhaloes at the high 
mass end, with Mud > 10 9 Mg, of all simulations (colour 
coded). The vertical dashed line marks the 25 kpc radius 
and all the objects to the left of this line are our crossed ob¬ 
jects. Instead of using a simple distance cut (and a mass cut 
at M t id > 10 8 Mg), the strength of the tidal forces on the 
disc and resonances of the orbital motion and the disc could 
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100 


r peri [kpc] 


Figure 3. M t ; d against pericentre distance r per ; for all 7 simula¬ 
tions (colour coded). The vertical dashed line marks the 25 kpc 
crossing criteria radii; all the subhaloes to the left of this line are 
regarded as crossed candidates. The 4 dotted lines correspond to 
lines of same tidal force at 10 kpc from the galactic centre as 
objects with M t ; d of 1, 3, 10 and 100 X 10 s Mg (from right to 
left) at r per ;=25 kpc. All subhaloes to the left of these lines have 
higher tidal forces. 


tion in all seven simulations together with statistics of the 
“crossed” sub-sample. N su b is the total number of subhaloes 
in each simulation suite with M t id ^ 10 6 M@, while p au b 
shows the fraction of subhaloes within r 50 . The percentage 
fraction which cross the introduced 25 kpc sphere is repre¬ 
sented by Pcross- VLII and Aq-F2 have the largest p su b values 
together with the highest crossed fractions. However, Aq-E2 
possesses the least number of subhaloes that cross the disc, 
although it has the fourth highest substructure fraction in¬ 
side ?' 5 o. The last four columns of tablc[3]show the number 
of crossed subhaloes with masses ^ 1, 3, 5 and 10 x 10 8 
Mq. We found a range of subhaloes from 12 to 24 with M t id 
^ 1 x 10 8 Mq which provides a good sample for analysing 
the impact of both the number and mass range of satellites. 
Aq-D2 has the largest number of subhaloes with mass ^ 1 
x 10 9 Mq and Aq-F2 is the only simulation that has one 
candidate ) 1 x 10 10 Mq. 

The time evolution of the satellites’ orbits from our N- 
body simulations are shown in Fig. [4] for all 7 simulations. 
The colours correspond to the mass range of the satellites; 
10 8 ^ M t id < 5 x 10 8 (grey) and A/ t id > 5 x 1O 8 M 0 fol¬ 
lowing a colour map. The dashed lines represent the 25 kpc 
crossing radius. The real pericentric distances from our sim¬ 
ulations are very close to the estimated pericentres found 
in the analysis that used an orbit integrator and the back¬ 
ground potential of the host halo (section [2.2| ). The presence 
of additional mass from bulge+disc slightly reduces r per i, 
hence the effective limit of selected satellites is < 25 kpc. 
The crossing of satellites is rather continuous with time. 


have been used. The four dotted lines correspond to lines 
of constant tidal force at 10 kpc from the galactic centre. 
These correspond to objects with M t id of 1, 3, 10 and 100 x 
10 8 Mq (from right to left) at r per i=25kpc. All sublialoes to 
the left of these lines have higher tidal forces. There is one 
object in Aq-E2 with a mass of ~ 1 O 1 O M 0 and a pericen¬ 
tre at 28 kpc, to the left of the red line (10 9 Mq) which has 
a tidal force greater than a subhalo with M t id=lO 9 M 0 at 
r pe ri=25 kpc and should have been included; however there 
are 2 subhaloes in Aq-E2 with larger tidal force (left of green 
line). It turned out, that only satellites left of the green line 
contribute to the disc heating on a measurable level. We 
conclude that we did not miss a significant fraction of inter¬ 
actions, which would change the derived dynamical heating 
of the disc. 

Previous analytical studies proposed a linear scaling be¬ 
tween the heating of the galactic disc due to satellites and 
the subhalo mass M aa t- However |Hopkins et al.| ( |2008| ) in¬ 
vestigated the importance of such dependence using "live” 
components for disc and subhaloes together with realistic 
non-circular orbits of satellites. They concluded that the 
heating rate per Gyr due to continuous infall, which forms 
the basis of our study, scales non-linearly as oc M 2 at tak¬ 
ing into account dynamical friction. This result contradicts 
earlier work which had only focused on the single interac¬ 
tion of subhaloes with discs and scaled linearly as oc M aa t 
per merging event. In other words, most of the impact is ex¬ 
pected from the most massive subhaloes found in the scatter 
in the high-mass end of the subhalo mass spectrum (Fig. [ 2 ]). 

I 11 this study we focus on the interaction of subhaloes 
with Mtid ^ 10 8 Mq corresponding to M t id/ M d ; ac > 0.003. 
Table|3] includes information regarding the subhalo distribu- 


2.3 Subhalo statistics 

Prior to the numerical analysis of the interaction of satellites 
and the galactic thin disc it is useful to discuss the properties 
of all subhaloes in our seven cosmological simulation suites. 
Fig-0 represents the inverse cumulative histogram of three 
quantities: N(>M) (dashed), M t id (solid) and M 2 d (dash 
dotted) for all the crossed subhaloes initially within rso. 
Aquarius simulations have a similar range for N au b. Aq-F2 
shows a significant jump at the high mass end of M t id due 
to the presence of one subhalo with M t id = 5.9 x 10 10 Mq 
which actually resides in the minor merger range. For Aq- 
F2 80% of the mass resides above 10 9 M 0 , while for the 
rest of the simulations the same percentage is reached at 
10' Mq. As mentioned above the expected impact strength 
scales with M t 2 id ; more than 95% of the impact is expected 
from the really massive subhalo for Aq-F2. Aq-D2 has the 
second highest expected impact due to massive subhaloes 
and Aq-C2 has the lowest. In all cases more than 90% of the 
impact is expected from subhaloes above 1 O 8 M 0 . 

The cumulative histogram of V ma x for all (solid) and 
crossed (dashed) subhaloes is shown in Fig.[6j The V max val¬ 
ues are normalised using the circular velocity (V 50 ) of the 
host halo at rso. About 90% of the subhaloes have V max < 
O.O 5 V 50 since the simulations are dominated by lower mass 
subhaloes. Moving towards subhaloes with V max > O.IV 50 , 
the scatter between the simulations increases due to low 
number statistics. For all the subhaloes, VLII lies on aver¬ 
age below Aquarius over the whole spectrum. For the crossed 
subhaloes it appears the deviation between VLII and Aquar¬ 
ius becomes less significant. The two dashed lines correspond 
to the average of the best-fitting power laws over all the 
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Table 3. Number statistics of subhaloes; p su b is the percentage of the subhaloes originally within r $o while p c ross shows the percentage 
of subhaloes inside r §q which cross the disc in 2 Gyr. ★: One of the subhaloes has M > 10 10 Mq 


Simulation 

N su b 

Psub 

[%] 

Pcross 

[%] 

Ncross 

(> 10 8 Mq) 

N 

■‘■’•cross 

(> 3 X 10 8 Mq) 

Ncross 

(> 5 X 10 8 Mq) 

Ncross 

(> 10 9 Mq) 

VLII 

28,618 

30.36 

12.71 

21 

5 

3 

1 
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Figure 5. Inverse cumulative histogram of number N (dashed), 
M t ^ (solid) and M 2 id (dotted) of crossed subhaloes for all our 
seven simulations as function of subhalo mass M t The vertical 
line marks the mass cut at 10 6 Mq. 


seven simulations except in the case of all the subhaloes 
were we have excluded VLII. The fitted power laws have the 
form 

N(> V max ) - P (Vmax/VsO)“ ( 6 ) 


/3 = 0.176 x 10 s a = -2.96 all 

/3 = 0.035 x 10 4 a = -2.87 crossed. 


The slope of the power laws are similar, indicating that the 
crossed subhaloes are representative of the total population. 
Also the slope agrees well with that mentioned in |Springel| 
et al. ( 2008| l. 


Fig.[T| presents two quantities p C ross,i (top) and p C rosa,2 
(bottom) which are defined as N cross (? - )/N a n(r) and 
N crosa (r)/N cross , respectively. The crossed subhaloes located 
at r > 100 kpc contribute to less than 20% of total popula¬ 
tion at such distance. However these subhaloes make up ~ 
50% of the total crossed fraction p C ross, 2 . 

For the purpose of numerical analysis each crossed sub¬ 
halo candidate will be generated as a distribution of particles 
following a NFW profile. To reproduce these subhaloes with 
a corresponding NFW density profile we determine the ra¬ 
dius, r t id, enclosing the provided mass Mad, using the scale 



Figure 6. Cumulative Vmax/Vso histogram for all (solid) and 
crossed (dashed) subhaloes; here V m ax is the maximum circular 
velocity of the subhalo while, V 50 is the circular velocity of the 
main halo at r 50 . 


radius r s and V max of subhaloes. Using 

TVmax = 2.163 r s , (7) 


/ Vmax \ —4 iG/) sat (r v „ x ) (8) 

\ f 'vmax / 

and 


Mad 


4 7t p s rl 


ln( 1 + ctid) 


Ctid 

1 + Ctid 


(9) 


together with equation 0 outline the most important steps 
in determining the tidal radius r t id of subhaloes via calcu¬ 
lating ctid which is the ratio of rad/r s . The median of the 
concentration c t id as a function of distance from the centre 
of the host halo is shown in Fig. [8] for all (top) and crossed 
(bottom) subhaloes. The outer most bin at r ~ rso ex¬ 
tends to infinity. For guiding the eye, best power law fits 
for all subhaloes (dashed) and for the crossed ones (dot¬ 
ted) are shown. In both populations we observe higher ctid 
with increasing distance from the host. Hence the subhaloes 
located further away from the centre tend towards higher 
concentrations. This trend is due to the tidal mass loss of 
the subhaloes and is well known ( fNagai fc Kravt sov 2005). 
The scatter between different simulations is insignificant. It 
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Figure 4. The time evolution of the orbits of subhaloes from our N-body simulations in spherical coordinates for all seven host haloes. 

Subhaloes with 10 8 . M t i<j < 5 X 1O 8 M0 are shown using grey colour while objects with A-ftid ^ 5 X 10 8 Afg. are represented following 

a colour map. 


is notable that the concentration of subhaloes located at r 
> 30kpc are, on average, lower for the crossed sample as 
seen by the deviation of the slope between the fitted power 
laws. At r=150kpc where the maximum p c roS s ,2 is observed 
(Fig.[7|, the total population has higher Ctid by a factor of 
1.12. For r < 30kpc the two populations show very com¬ 
parable concentrations. This suggests that subhaloes on ra¬ 
dial orbits loose more mass, when entering the main halo 
at rso already. Fig.[9]is the normalised histogram of c t id for 
all (solid) and crossed (dashed) subhaloes. The peak of the 


distribution for the total population is around Ctid ~ 4.5 
while for the crossed subhaloes we observe a shift of the 
distribution to lower concentrations with the peak around 
3.6. This indicates that lower concentration is favoured by 
crossed satellites. The reason for this behaviour consists of 
two points. The mean concentration of crossed subhaloes is 
smaller for r > 100 kpc (Fig.[8]l and the fraction of crossed 
subhaloes, p C ross,i, is decreasing with increasing r leading to 
a biased selection towards lower cud- 

The median of the concentration as function of sub- 
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log(r [kpc]) 

Figure 7. Top : Fraction of crossed to all subhaloes p C ro BB ,l 
for spherical shells from centre of their host halo, defined as 
N C ro SS (r)/N all (r). Bottom : Normalised radial distribution of all 
crossed subhaloes, p croBB , 2 , representing N crosB (r)/N croBB . 


halo mass for all (top) and crossed (bottom) subhaloes is 
plotted in Fig. |10| Here and also for Fig.[l3]the final bin in¬ 
cludes all the subhaloes with larger mass. The scatter for the 
crossed subhaloes is smaller than for the total population. 
The dashed and dotted lines represent best fits for all and 
crossed samples following a power law 


Ctid — C8 


( Mtid 
= CR — 


Vio s M 0 y ’ 

where cs corresponds to the Ctid value at tidal mass of 
10 8 Mq. The slope of the power law is identical in both 
plots (77 = 0.05). We observe a shift to lower concentrations 
for crossed subhaloes by a factor of 1.35 which corresponds 
to the shifted maximum in Fig. [9] The agreement between 
VLII and Aquarius is particularly noticeable for the crossed 
populations. 

represents the distribution of pericentres (solid) 


Fig- 


11 


and initial distances (dashed) of subhaloes for all seven sim¬ 
ulations. The vertical lines mark the position of the 25 kpc 
crossing criterion. VLII has the lowest number of subhaloes 
both with pericentres < 25 kpc and also over the whole dis¬ 
tance range. This is partially due to the fact that VLII is 
dominated by low-mass subhaloes, with M t id ^ 10 6 Mq , 
which are ignored in this analysis. There exists some varia¬ 
tion in the maximum of the pericentre distribution between 



!og(r [kpc]) 

Figure 8. Mean concentration Ctid against the distance from 
centre of the host halo, r, for all (top) and crossed (bottom) sub¬ 
haloes. The dashed and dotted lines represent the average best¬ 
fitting power laws for the total and crossed candidates, respec¬ 
tively. 


the simulations, but most subhalo pericentres lie in the range 
10 < r < 100 kpc. 

Another interesting aspect is the angular distribution of 
subhaloes in the simulations. The cosine of the inclination 
angle i (solid) and the azimuthal angle <j> (dashed) of initial 
positions are represented in Fig. [12] which are measures of 
isotropy. These quantities are calculated by 


cos (i) 


2 

r 


and 


( 11 ) 


</> = tan -1 (!) - 180 < <f> < 180, (12) 

where x,y and z represent the position of the subhaloes 
in Cartesian coordinates and r is the distance to the host 
halo’s centre in spherical coordinates. 

It is important to note in this analysis, the x,y and 
z directions do not necessarily imply any favoured orienta¬ 
tion of the system but rather the chosen disc orientation 
(section[3|. In the case of cos (i), the inclination angle i is 
measured with respect to the 2 -component. Therefore cos(*) 
of 1 corresponds to the object which is located on the posi¬ 
tive 2 -axis and vice versa for cos(i) of -1. For the Aquarius 
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Figure 9. Normalised distribution of Ctid for all (solid) and 
crossed (dashed) subhaloes. 



log(M tia [M 0 ]) 


Figure 10. Concentration c t id against f° r all (top) and 

crossed (bottom) subhaloes. The dashed and dotted lines repre¬ 
sents the average best-fitting power laws for the total and crossed 
candidates. 



r [kpc] 


Figure 11. Pericentre (solid) and initial distance (dashed) his¬ 
tograms for all subhaloes. The vertical dotted lines mark the 
25 kpc crossing radii. 


simulations this quantity has a rather homogeneous distri¬ 
bution with Aq-E2 showing a slight peak around 0, indicat¬ 
ing that more subhaloes are concentrated towards the x-y 
plane. However, it is interesting that in the case of VLII 
we clearly see distinct peaks close to 1 and -1, suggesting 
that most of the subhaloes are close to the 2 -axis, i.e. on 
more polar orbits with respect to the disc. The azimuthal 
angle </> tells us how the sublialoes are distributed initially 
when observed from a top-down view (on the x-y plane). We 
observe the largest deviation from an isotropic distribution 
for VLII and Aq-A2. The bimodal behaviour indicates that 
more subhaloes are located around ± 90° corresponding to 
the y-axis. Aq-B2 has the most isotropic distribution. We 
also performed a similar analysis for the crossed subhaloes 
and concluded that for the Aquarius simulations the distri¬ 
butions for both quantities become more isotropic - flat in 
cos(i) and 4> - while for VLII the bimodality is still found. 

One of the important quantities in order to characterize 
the shape of orbits is the eccentricity e which in calculated 
using 


e— Jl —— where fi = G(M t - ld + M h (r)), (13) 

V 
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Figure 12. Distribution in cos(i) (solid) and <j> (dashed) for all 
subhaloes; where i is the subhalo inclination angle with respect 
to the ^-component and <j> is the azimuthal angle. 


I = I r X v I. 


(14) 


a = -and 

2e 


(15) 


e = < 16 > 

In the above equation, fi corresponds to the sum of the stan¬ 
dard gravitational parameters of the subhalo and the host 
halo with enclosed mass Mh at distance r. The specific an¬ 
gular momentum l is calculated using the velocity and posi¬ 
tion of the subhalo with respect to the host. The semi-major 
axis a and the specific orbital energy e are related through 


equation (161. Only in a Kepler potential, the eccentricity 


and semi-major axis are well-defined and constant along the 
orbit. A general potential using the enclosed mass, gives 
an upper limit for the eccentricities and underestimates the 
pericentre distances. 

Fig. [13] shows the median eccentricity against M t id for 
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Figure 13. Eccentricity e against M,;,j for all (top) and crossed 
(bottom) subhaloes. 


all (top) and crossed (bottom) subhaloes. VLII possesses 
slightly higher e values than Aquarius over the whole mass 
range for the total population. Aquarius simulations have 
a relatively similar behaviour in the M t id < 10 9 Mq zone, 
while the scatter in the higher mass regime is expected. For 
the crossed subhaloes, the difference between the simulations 
is minimized and there exists a minor increase of e with 
Mtid- In other words the orbits of crossed subhaloes become 
more parabolic, closer to e = 1, for higher masses. 

More than 85% of the subhaloes in all seven simulations 
have 0.5 < e < 1 - these subhaloes have bound elliptic orbits 
with respect to the enclosed halo mass at the corresponding 
radii. Also ^ 10% have hyperbolic trajectories with e > 1 
and regarded as candidates which would cross once as they 
pass by the galactic disc. VLII has the highest fraction of 
subhaloes with unbound orbits. 

In order to determine how homogeneous is the infall 
time distribution of subhaloes, we looked at the pericentric 
passage time of the subhaloes which correspond to the time 
of the closest approach to the centre of the host halo. The ex¬ 
pected pericentric time of crossed subhaloes for all the seven 
simulations is shown in Fig |l 4| The distribution is fairly ho¬ 
mogeneous while we observe a shallow peak around 1 Gyr 
time; hence the impact of the subhaloes is expected to be 
continuous with time. There is a lack of encounters for the 
first 0.5 Gyr which might be due to mass loss of subhaloes al- 





























^■pericentre [MyI"] 

Figure 14. Distribution of pericentric time t per i centre , time of 
reaching the closest pericentre, for crossed subhaloes. Multiple 
pericentre passages are not included. 


ready located at close distances. The possible impact of this 
delay on heating of the galactic disc is discussed in section]!] 
Also the drop seen in all the simulations after 1.5 Gyr could 
be the result of simulation box’s finite volume size. This 
plot can robustly tell us about the interaction rate per Gyr 
of subhaloes with the host. 

The VLII and Aquarius simulations do not show sig¬ 
nificant differences in their spatial distribution. The higher 
mean eccentricities of the VLII satellites is the main differ¬ 
ence compared to the Aquarius simulations. However for the 
crossed subhaloes such a difference disappears. 


3 GALAXY N-BODY MODEL 

Prior to the numerical analysis we need to set up our desired 
galaxy models as a distribution of particles which include 
disc, bulge, DM halo and also the satellite component. As 
mentioned earlier, one of the most important aims of this 
work is to resolve the dynamical properties of the galactic 
disc in the absence of any bias originating from a system 
which is not originally in equilibrium. 


3.1 Multi-component Milky Way galaxy 


In the past two decades there has been a number of at¬ 
tempts to generate initial conditions for multi-component 
spiral galaxies that include a disc, bulge and DM halo com¬ 
ponent, and producing a distribution of particles which fol¬ 
low particular density profiles. We have tested a few such 
attempts, however the true equilibrium and steady state of 
the disc component is questionable. For this study we take 
advantage of a recent code (GallC) which employs a different 
approach to the made-to-measure method: constructing ICs 
for multi-component axisymmetric systems via an iterative 
scheme (Yurin & Springel [2014 1. 

The stellar disc component has an exponential radial 
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profile while for the vertical direction an isothermal profile 
is adopted 


/ r* \ Md 

Pd(i? - )= 4W? eXP 


sech 2 ( — 

h I \ z 0 


(17) 


In the above equation h and zo correspond to the radial 
scale length and the thickness of the disc; the total mass of 
the disc is represented by M d ■ For the scaling of the disc 
surface density, we use the value of the thin disc at the 
solar neighbourhood E so i proposed by Just & Jahreifi (2010) 
which is then converted into Md via 


M d = 2irY, aol h 2 exp ( ^ 


(18) 


where Ro = 8 kpc is supported by Gillessen et al. (2009 ) 
and |Reid et al.| |2014| ). The value for the radial scale length 
h is chosen to be 2.8 kpc. Since the disc mass is not sensitive 
to the scale length over a large range, when fixing the local 
surface density, a strong influence from smaller scale lengths 
as discussed in recent literature is not expected. The GallC 
code allows us to constrain the kinematical preferences of 
axisymmetric components. For the disc, we use the ratio of 
the radial and vertical velocity dispersions (Jr/(j z which con¬ 
verges to a unique solution during the IC construction. The 
bulge follows a spherically symmetric and isotropic Hern- 
quist profile (equation |20| ) . The DM halo is set up in two 
stages. The first step consists of using properties of our DM 
halo such as A, c and V 200 (the spin parameter, concen¬ 
tration and total circular velocity of the system at )" 2 oo) 
extracted from the seven simulations. These quantities to¬ 
gether with the employed cosmology allows us to infer the 
total mass (disc+bulge+DM halo) of the system, M to t, using 


Mtot — 


10 GHo 


(19) 


The spin parameter A, which determines the size of the disc, 
possesses a range of values for different haloes and is depen¬ 
dent on the concentration. Since we would like to have sim¬ 
ilar disc profiles among all simulations, this quantity is cal¬ 
culated such that the discs have similar scale lengths. These 
systems lie in the low spin parameter range. The code takes 
a reverse procedure and using the inputted A together with 
the disc mass fraction m d, it calculates the scale length of 
the disc (e.g. IMo, Mao fe White||l998 Springel, Di Matteo 


& Hernquist 20051. In the second step, the code uses the 


technique where a matching Hernquist profile is created for 
our desired NFW halo. This profile has a similar density 
profile to NFW in the inner part while falling off steeper in 
the outer regions according to 

M hern a 


PheTnlX ) — 


2-7T r(a + r) 3 ’ 


( 20 ) 


where Mhern is the total mass of the halo and a is the scale 
radius of the Hernquist profile. |Springel, Di Matteo fe Hern-| 
quist (2005) introduced a method for constructing a Hern¬ 


quist model with a density profile which matches that of a 
particular NFW profile. This is achieved by associating a 
Hernquist DM halo with total mass - Mhern - same as the 
NFW halo mass within r 2 oo where scale lengths are related 
via the NFW halo concentration 

a = r B y/2[ln(l + c) — c/(l + c)]. 


( 21 ) 
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Table 4. The characteristics of the components for the Milky 
Way. N represents the number of particles for the component 
while the subscripts d,b and h correspond to the disc, bulge and 
the halo. The solar neighbourhood surface density £ ao i together 
with the thin disc scale height and scale radius, zq and h corre¬ 
spond to the model proposed by|Just &; JahreTE] (|20 1 0~|. 


Quantity 

Value 

Ssol 

39.9 Mqpc ~' 2 

h 

2.8 kpc 

zo 

277 pc 

M d 

3.42 x 10 10 Mq 

Cfi/Cz 

1.5 

N d 

10 million 

a b 

0.35 kpc 

M b 

1.9 X 10 10 Mq 

N 6 

0.5 million 

A?200 

1.77 x 10 12 M @ 

A 

0.021-0.028 

N h 

4 million 


We used this approach to generate the matching Hernquist 
potential for the parent haloes and employ them in the or¬ 
bit integration analysis. We also found negligible deviations 
between the orbits of subhaloes in the NFW and the cor¬ 
responding Hernquist haloes of the order of < few percent. 
All the DM haloes in this study are generated with a spher¬ 
ical density profile for simplicity, while the original profile 
from the cosmological simulations are non-spherical, i.e. pro- 
late/oblate/triaxial, (Allgood et al.|2006 l. 

Some of the quantities used for the set up of the initial 
conditions are listed in tablcjd] The A ranges from 0.021- 
0.028 depending on the simulation and the DM halo which 
is constructed. Each system contains 10, 0.5 and 4 million 
particles for the disc, bulge and the halo respectively. Ac¬ 
cording to Khoperskov et al. (20071 an insufficient number 
of disc particles could damp out the modes by introducing 
noise; hence 10 million disc particles allows for a physically 
reliable analysis of the disc. 


3.2 Satellite model 

The satellites with M t id ^ 10 s Mq are generated using their 
Myd, c t id and r s calculated in section [273| tog ether with the 
distribution function introduced by Widrow (20001 and im¬ 
plemented in Lora et al. (20131 for a cuspy NFW profile. 
We need to employ a cut-off radius for the satellites due 
to the infinite cumulative mass of NFW profile at r —► oo; 
hence rtid was chosen. All the satellites for this study are 
represented using 50,000 particles. 


3.3 N-body code 


The numerical simulations in this paper are performed 
using the efficient parallelised N-body particle-mesh code 
SUPERBOX-10 (Bien, Brandt & .J ust| [2013). Since the in¬ 
dividual gravitational interaction of particles from simi¬ 
lar/different components, also known as two-body relax¬ 
ation, has a time-scale of the order of the Hubble time we 
can regard them as collisionless ( Binney fe Tremaine|2008 1. 
SUPERBOX-10 has three nested grids at different resolution 
for each galaxy, where the inner most grid with the highest 
resolution allows the inner dense structure of every galaxy 


to be resolved and the outer most grid contains the whole 
simulation and coincides for all galaxies. Each grid is divided 
into cubic cells and we decided to use n 3 cells with n=256 
instead of 128 in order to reduce the energy conservation 
error to < 0.4% after 4 Gyr of evolution. The size of each 
cell depends on the radius of the grid R gr id,i and calculated 
as 


dcell 


2 Rgrid,i 

n — 4 


( 22 ) 


The outermost grid contains all the particles in the simu¬ 
lation and hence has a size of few hundred kpc. The cal¬ 
culation of the potential is performed via transforming the 
cell densities into potential using parallelised Fast Fourier 
Transform (FFT) together with Green’s function. The com¬ 
putation time scales as NlogN in the case of FFT while 
direct calculation scales as N 2 . In addition, the orbital in¬ 
tegration of the particles are done using a leapfrog method 
and a constant time step (Fellhauer et al. 2000). It is im¬ 
portant to choose an appropriate time step as it plays an 
essential role in the force calculations; it needs to be short 
enough that a typical particle in the inner most grid moves 
a maximum of one cell during the step, which was chosen to 
be 0.1 Myr. 

In order to enhance the resolving power of the vertical 
structure of the galactic disc, SUPERBOX-10 takes advan¬ 
tage of a flattening parameter q which flattens the middle 
grid along the 2 -axis. This technique improves the calcu¬ 
lation of the acceleration in the vertical direction. For all 
our simulations we use similar grid radii for the Milky Way 
galaxy; Ri nn = 3.5kpc, R mid = 30kpc and R out =500kpc. 
The cell size of each grid are di nn = 0.028 kpc, d m id = 0.24 kpc 
and d ou t = 4.0kpc. The flattening factor q has the value of 
0.3, which increases the resolution of the middle grid in the 2 
direction by a factor of ~ 3.3. According to the work by|Just| 


& Khan 2011) the effective resolution in the force calcula¬ 


tions using the particle-mesh code is d e fi = d ce ii/2. Hence 
we reach a resolution of 40 parsecs in the vertical direction 
of the entire disc in the middle grid. 

Each satellite has its own 3 grids, where the sizes of the 
inner and middle grids are chosen as 2? - s and 2r t id; the outer 
grid has the same size as the MW system, 500 kpc. 


3.4 Analysis of the disc 

The disc is subject to being shifted, tilted and distorted with 
respect to the original coordinate system due to evolution 
and mostly the bombardment by the satellites. Therefore we 
defined at each time-step a new coordinate system for the 
analysis of the disc structure by moving the origin to the 
density centre of the host galaxy. The velocities of the disc 
particles are corrected for the density centre velocity, then 
the coordinate system is rotated to correct for the tilt. It 
turned out that the angular momentum vector of the disc 
is the most robust measure of the disc rotation axis. The 
positive 2 -axis is rotated to the angular momentum vector 
of the disc, taking into account all disc particles. Since we are 
interested in the radial range of 4-15 kpc (for an exponential 
disc with flat rotation curve 95% of the disc mass and 87% 
of the angular momentum is enclosed within 5 h= 14kpc). 
The tilt angle ranges from <1° for most of our simulations 
to ~ 4° in the case of Aq-F2 with satellites. 
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For the analysis of the disc structure we study the 
temporal evolution of the solar neighbourhood defined by 
7.5<R<8.5 kpc, and the radial profile at the end of the sim¬ 
ulation. We evaluate in radial rings, z me an which provides a 
measure of the mean value of the disc z-coordinate, in order 
to check the disc orientation and identify bending modes. 
We have checked that the inner disc is well-aligned with the 
new coordinate system in all cases. The root-mean-square of 
the disc z-coordinate in the solar neighbourhood calculated 
via 

Zrms — \/ + (z Zmean ) 2 +, |z Zmean | + 2 kpC (23) 

which is a measure of the disc thickening and the square of 
the vertical velocity dispersion are measures of disc’s vertical 
heating. We neglect any disc particle with |z — z me an| ^ 2 kpc 
for this analysis in order to avoid any outliers. 


3.5 Isolated system 


The multi-component Milky Way galaxy was evolved in iso¬ 
lation - in the absence of satellites - for 2+2 Gyr in the 
case of all 7 simulations. The initial 2 Gyr is for the purpose 
of reaching an equilibrium state while the latter 2 Gyr is 
to have a control group for comparison to the systems with 
satellites. The isolated evolution allows the components to 
adapt to their environment and also reach the equilibrium 
and steady state. Fig.[l5] shows the time evolution of 10- 
90% Lagrange radii of the Aq-B2 system (disc+bulge+halo) 
in isolation. This quantity represents the spherical radii at 
which 10-90% of the total mass are enclosed. The inner 50% 
of the system that extends beyond the disc does not show 
signs of expansion/contraction during 4 Gyr. The outer re¬ 
gions of the halo have only < 6% of expansion, which is 
expected due to the cut-off as the outer parts are loosely 
bound and few particles could in principle escape. 

The number density plot of the disc from the Aq-B2 
simulation is represented in Fig.[l6] where starting at top- 
left we observe the system at -2, 0, 1 and 2 Gyr. After 4 
Gyr of evolution the disc appears in steady state and we do 
not observe any strong instabilities within the disc. Also the 
red crosses correspond to centre of mass of the bulge which 
coincides with the disc centre of mass. 

Fig. 1 1 7| shows a few quantities describing the dynamics 
of the isolated system Aq-B2 - a fair representative of our 
simulation sample. In the top two rows different line colours 
correspond to the system at -2, -1, 0, 1 and 2 Gyr of evolu¬ 
tion as function of radial distance across the disc. Starting 
from the top left, the surface density E of the disc is repre¬ 
sented and stays unchanged throughout 4 Gyr of evolution 
in isolation. We observe some fluctuations at the outer disc, 
R > 10 kpc, in the initial 1 Gyr which dampens afterwards. 
The next plot shows the Toomre parameter Q, a measure of 
the disc stability and is calculated as 


ORK 

“ 3.36GE ’ 


(24) 


where <jr and k represent the radial velocity dispersion and 
the epicyc.lic frequency of the disc particles. Initially the disc 
has a Q value of 1.5 at the solar neighbourhood increasing 
to ~ 2.5 after few hundred Myr. The initial increase of or 
is still one of the caveats of the GallC code. The vertical 
velocity dispersion a, shows an increase of 45% at R = 15 kpc 



time [Gyr] 


Figure 15. 10-90% Lagrange radii of the isolated MW for Aq-B2 
simulation. The 10% radii corresponds to the lowest value (dark 
purple) while 90% has the highest value (red) and the remaining 
colours in between show the 20-80% Lagrange radii from bottom 
to top. 


after 4 Gyr while for the solar neighbourhood we observe 
an increase of < 12% only. Accordingly, z mean increases by 
maximum of 25 pc for the inner 8 kpc while the outer part 
experiences an increase of 100 pc after 4 Gyr. The next plot 
shows the time evolution of z rms in the solar neighbourhood. 
The increase of the disc thickness after 3 Gyr is ~ 15 % 
which reaches 320 pc while in the case of the cr 2 in the solar 
neighbourhood a vertical heating of 40km 2 s -2 - 16 % - is 
seen in the last 3 Gyr. 


3.6 Fiill simulations 

After obtaining an isolated system which satisfies the equi¬ 
librium criteria, we add the subhalo candidates using their 
extracted initial positions and velocities at t=0 Gyr. Since 
there is no massive satellite with M t id > 5 x 10 8 A7 q closer 
than 20 kpc (see Fig. [4]), we don’t expect an unphysical reac¬ 
tion of the disc by the sudden implementation of the satel¬ 
lites. Each of the seven systems were simulated for 2 Gyr. 


4 RESULTS 

In this section we discuss the results obtained for the analysis 
of satellites’ impact on the disc. We focus on the vertical 
heating of the galactic disc due to the infalling satellites, 
expressed in terms of thickening the disc and increasing the 
vertical velocity dispersion of disc particles. In the first part 
we discuss the heating in the solar neighbourhood, 7.5 kpc 
< R < 8.5 kpc, for the duration of 2 Gyr while in the second 
part we analyse the impact of satellites across the whole disc 
radial range. In the last part, we argue about the relevance 
of the orientation of satellites’ distribution to the vertical 
heating measurements for two cases of Aq-D2 rotated along 
the x and y axes by 90 degrees. 
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Figure 16. The number density plot for the isolated system in x-y, x-z and y-z view. From top left: Aq-B2 simulation at initial time (-2 
Gyr), 0, 1 and 2 Gyr of evolution. The red cross represents the centre of the mass of the bulge component. 


4.1 Solar neighbourhood 


In this study, with advantage of our high resolution analysis 
of the Milky Way disc, we are able to resolve the structure 
of the disc at 40 pc scale together with the large number 
of disc particles. The next results focus on the time evolu¬ 
tion of the 2 rma and vertical velocity dispersion in the solar 
neighbourhood for all the particles lying in the radial range 
7.5 - 8.5 kpc from the disc centre. The top panel of Fig. |18| 
shows the time evolution of the root-mean-square of the disc 
height subtracted by the initial values, z r ms,t — 2 rms ,o for sys¬ 
tems with satellites (solid) and isolated (dashed) at the solar 
neighbourhood for VLII (blue), Aq-B2 (green), Aq-D2 (red) 
and Aq-F2 (orange). During the initial 0.4 Gyr the thick¬ 
ening is minimal since no massive subhalo with M t id ^ 10 9 
M@ crosses the disc within this period (Fig.[4j). The reason 
for this lack of interaction in the first 0.4 Gyr may be the 
result of initial conditions at 2 = 0 , since we did not take into 
account the tidal mass loss which already happened for satel¬ 


lites in the inner region of the host halo (Nagai & Kravtsov 
2005). In order to take this into account; we will compare 
the dynamical evolution in the simulations in addition to 
a delayed heating by 0.5 Gyr. Aq-F2 reaches a maximum 
thickness difference of 75 pc (~ twice as much as in the iso¬ 
lated case at 2 Gyr); while in the remaining simulations the 
thickness increases by 30 pc at most. 

As mentioned in section[l] the primary aim of this study 
is to understand the contribution of realistic cosmological 
satellites to the heating of the galactic disc. This process can 
be best quantified using the squared vertical dispersion cr? of 
the disc particles. The most reliable measurement of the ver¬ 
tical heating of the Milky Way disc has been obtained for the 
solar neighbourhood (e.g. Hol mberg et al.|20 07 2009) with 
the value of d(cr^)/dt = 72knT : s - ^ Gyr -1 following a power 
law with exponent 27 = 1.06 for a\ which means a constant 
heating rate. The time evolution of the quantity t - <r ^ 0 
(the difference between the system at time t and the initial 
time when the satellites are added) for systems with satel- 
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Figure 17. The properties of the disc in the isolated Aq-B2 simulation. From the top left as function of radial distance from the 
centre of disc: Surface density of the disc E, Toomre parameter Q, vertical velocity dispersion a z , vertical cri ,, mean disc height z m ean, 
root-mean-square thickness 2 rms . The lines represent the results at -2, -1, 0, 1 and 2 Gyr. The two bottom figures show the 2 rms and a z 
in the solar neighbourhood, 7.5 < R < 8.5, as function of time. 


lites and isolated in the sola r neighbourhood is presented in 
the bottom panel of Fig. 18 with <t? j0 = 260 km 2 s -2 . As with 
the evolution of z rms we observe a more significant heating 
in the case of Aq-F2 - 25% at 2 Gyr compared to initial 
while only very small fluctuations are seen for the remain¬ 
ing simulations. The self-heating of the disc is very small, 
on scale of ~ 10 km 2 s -2 Gyr -1 . The dotted and dot-dashed 
lines, respectively, represent the observed vertical heating 
rate in the solar neighbourhood and the observed value al¬ 
lowing for 0.5 Gyr delay which accounts for possible initial 
underestimation of heating due to mass loss of subhaloes al¬ 
ready located in the inner halo. The jumps observed in the 
velocity dispersion due to massive mergers are seen in other 
simulations where such mergers cause visible signatures in 


the velocity dispersion profiles (Martig, Minchev & Flynn 


2013). 


Aq-D2 has twice as many subhaloes with Mtid ^ 10 9 
Mq compared to Aq-B2 which, in theory, translates into a 
larger impact. We also observe twice as much heating com¬ 
pared to the isolated cases between the two simulations in 
the solar neighbourhood. 


The mean values of 2 rms ,t — z rm s,o for all 7 simulations 
with satellites (solid), with satellites except Aq-F2 (dotted) 
and isolated system (dashed) are shown in the top panel 
of Fig. [T9] with their subsequent 1 <t statistical dispersion 
around the mean. The blue, pink and green colours cor¬ 
respond to all the simulations with satellites, with satellites 
except Aq-F2, and for the isolated system, respectively. We 
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Figure 18. Top: Time evolution of the root-mean-square of the 
disc thickness subtracted by the initial value for systems with 
satellites (solid) and isolated (dashed) for VLII (blue), Aq-B2 
(green), Aq-D2 (red) and Aq-F2 (orange) simulations in the so¬ 
lar neighbourhood. Bottom: The time evolution of the difference 
of the disc vertical heating for VLII (blue), Aq-B2 (green), Aq- 
D2 (red) and Aq-F2 (orange) simulations. The dotted and dot- 
dashed lines represent the observed vertical heating rate in the 
solar neighbourhood and the observed value with 0.5 Gyr delay, 
respectively. 


observe a very shallow increase of the thickness with time in 
all systems after the initial 0.4 Gyr. The mean value of all 
the systems with satellites experiences an increase of 35 pc 
compared to the initial state and only < 10 pc - ~ 40% - 
compared to the isolated system. If we exclude the impact 
from Aq-F2, the thickening in the solar neighbourhood from 
satellites is almost indistinguishable from the isolated case. 

According to the bottom panel of Fig. |19| there exists a 
slow increasing trend of the mean value of aT t - a^ 0 for all 
systems. The heating rate for the initial 1 Gyr are 20, 15 and 
10 km 2 s -2 Gyr -1 for systems with satellites, with satellites 
except Aq-F2, and the isolated system, respectively, which 
reduce to 15 km 2 s -2 Gyr -1 for all the systems later. This 
means the rate at which the vertical dispersion increases 


with time in the solar neighbourhood decreases slightly af¬ 
ter 1 Gyr. Not taking into account Aq-F2, the difference be¬ 
tween the mean heating due to satellites and the isolated 
systems is very small and only 5km 2 s -2 at 2 Gyr. The 
la dispersion region from all the simulations (blue) reaches 
50 kin 2 s -2 at the final snapshot which is 35% of the observed 
value of 144 km 2 s -2 . The additional heating for subhaloes 
below 10 10 Mq is of the order of ~ 10-15%. The dotted and 
dot-dashed lines represent the observed vertical heating rate 
in the solar neighbourhood and the observed value with 0.5 
Gyr delay, respectively. The Aq-F2 simulation suggests in 
order to reach the observed heating rate, that a few en¬ 
counters per Gyr of LMC type satellites with pericentres 
below 20 kpc are required. Our result fails by more than 3er 
to reproduce the observed heating rate. These results are 
consistent with the recent observational work by|Ruchti et| 


al. (2015). Using the Gaia-ESO spectroscopic survey, they 


found no evidence for a significant dark matter disc in the 
MW and that the Galactic disc has had a rather quiescent 
merger history for the past ~ 9 Gyr. 


4.2 The impact across the disc 

In addition to the solar neighbourhood, the global vertical 
heating across the whole disc was also investigated and is 
presented in this section. 

The root-mean-square of the disc height at 2 Gyr sub¬ 
tracted by the initial values, z rm s,t — Zrms,o ; as a function of 
radial distance from the disc centre for systems with satel¬ 
lites (solid) and isolated (dashed) for VLII (blue), Aq-B2 
(green), Aq-D2 (red) and Aq-F2 (orange) simulations are 
shown in the top panel of Fig. |20| this quantity is calcu¬ 
lated in radial bins of equal particle numbers. The vertical 
heating (oTt-cr^o) of the disc is shown in the bottom panel 
of Fig. |20| The solid and dashed lines illustrate the systems 
with satellites and the isolated control systems; the results 
correspond to the time of 2 Gyr which is the final snapshot 
of our simulations. The deviation between the control sys¬ 
tem and the system with satellites is mainly visible in the 
outer disc. In Aq-F2 the flare starts at the solar neighbour¬ 
hood and reaches a thickening of 550 pc at R= 15 kpc. For 
the VLII simulation we observe that the thickening of the 
isolated case is slightly larger than the case with satellites 
at some radii, however we still observe a positive increase 
of thickness compared to 0 Gyr time. As seen in previous 
studies |Ardi et al.] [2003 | we observe a small flaring in the 
outer parts of the disc which is amplified with the presence 
of the satellites; the particles which reside in the outer re¬ 
gions are less strongly bound to the disc and require less 
energy to be displaced. Aq-A2, Aq-C2 and Aq-E2 have the 
least thickening due to infalling satellites. 

The bottom panel of Fig. [20] illustrates the quantity 
Cz,t ~ a z, o- Aq-F2 dominates the heating - reaching ~ 
200 km 2 s -2 in the outer parts. The peak at R<4kpc cor¬ 
responds to a very small increase in velocity dispersion 
(from lOOkrns -1 to 102kms _1 at the centre). The ob¬ 
served vertical heating in the solar neighbourhood after 2 
Gyr (Holmberg et al. |2009[ ) is represented as a filled diamond 
while if we allow for 0.5 Gyr of delay we obtain the heat¬ 
ing represented by the empty diamond in Fig.[20]with val¬ 
ues of 144 and 108 km 2 s~ 2 , where 144 corresponds to twice 
the vertical heating rate of 72 km 2 s -2 Gyr -1 . The er 2 0 at 
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Figure 19. The mean difference of vertical heating between t = 2 
and t = 0 Gyr together with their la statistical dispersion region 
of all 7 systems with satellites (solid, blue), with satellites except 
Aq-F2 (dotted, pink) and isolated (dashed, green). The dotted 
and dot-dashed lines represent the observed vertical heating rate 
in the solar neighbourhood and the observed value with 0.5 Gyr 
delay. 




R [kpc] 

Figure 20. Top: The difference of root-mean-square of the disc 
height Zrms at 2 and 0 Gyr for systems with satellites (solid) and 
isolated (dashed) for VLII (blue), Aq-B2 (green), Aq-D2 (red) 
and Aq-F2 (orange). Bottom: The difference of the disc vertical 
heating between t = 2 and 0 Gyr for VLII (blue), Aq-B2 (green), 
Aq-D2 (red) and Aq-F2 (orange), filled and open diamonds rep¬ 
resent the observed value of vertical heating in the solar neigh¬ 
bourhood and the observed value with 0.5 Gyr delay. 


R = 15 kpc has a value of 40 km 2 s~ 2 and therefore the max¬ 
imum heating in the outer parts reaches ~ 500% compared 
to the initial value - excluding Aq-F2 we only have < 130%. 

Similar to the behaviour seen in the solar neighbour¬ 
hood, Aq-B2 and Aq-D2 possess comparable heating rates 
with 20 km 2 s -2 extra heating for Aq-D2. Aq-D2 has half of 
its pericentric passages during the initial 1.1 Gyr while the 
remaining passages occur during the final 0.1 Gyr. Hence it 
appears that the heating due to subhaloes with 10 9 < M t id 
M a < 10 10 Mq contributes only to ~ 40 km 2 s -2 of heating 
for R > 4 kpc. 

The top panel of Fig. 1 


21 


demonstrates the mean a 2 ;i 


- <j 2 j0 value of all seven simulations with satellites (solid), 
with satellites except Aq-F2 (dotted) and the isolated sys¬ 
tem (dashed), across the radial range at t = 2 Gyr. The 


shaded regions present the lcr statistical dispersion in each 
bin. The impact of Aq-F2 is seen as the difference of the 
dotted and solid lines which increases as we move towards 
outer regions. The mean disc height of all simulations with 
satellites increases by 150 pc after 2 Gyr at 15 kpc and in¬ 
cluding the dispersion this value reaches 340 pc; while the 
isolated systems only have an increase of 50 pc at this ra¬ 
dius, corresponding to self thickening of about 14%. Aq-F2 
increases the mean value by an additional 70 pc in the outer 
15 kpc regime. 

The mean value of the vertical heating between 0 and 
2 Gyr for all systems across the disc is shown in the bot¬ 
tom panel of Fig. |21| The shaded areas correspond to the 
lcr dispersion around the mean for the subsequent systems 
blue, pink and green. For all the systems with satellites, 
the mean value lies around 40 km 2 s~ 2 across the disc at 
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R>4kpc, while the simulations excluding Aq-F2 possess < 
20 km 2 s -2 in the outer parts. The solar neighbourhood heat¬ 
ing at R=8kpc of 144 km 2 s -2 and the value 108 km 2 s -2 
taking into account 0.5 Gyr delay are represented as filled 
and open diamonds. The lcr dispersion in the solar neigh¬ 
bourhood reaches only 50 km 2 s -2 . We observe a specific 
heating almost independent of galactocentric distances for 
R 4kpc. More energy is deposited to the inner parts; al¬ 
though, this region is not strongly perturbed due to both 
higher surface density and velocity dispersion. However, the 
outer parts with lower surface density experiences stronger 
relative heating together with flaring. 


4.3 Effect of satellites infall orientation 


In this section we present our results from the analysis of the 
impact of the orientation of satellites distribution relative to 
the disc on the measured vertical disc heating. For this pur¬ 
pose we only use the Aq-D2 simulation, which has the second 
most number of satellites candidates (23) above 10 s Mg and 
the most candidates with M t id ^ 10 9 Mg. This simulation 
is believed to be representative of our simulation suites and 


shows a typical level of anisotropy (Fig. 12 I. We decided to 


consider two additional cases of different orientation that 
consist of rotating the satellites’ frame by 90 degrees along 
the x axis, hereafter case B, and also performing the similar 
rotation along the y axis, case C. The original orientation is 
regarded as case A. These rotations correspond to changing 
the disc original axis of rotation, z, to rotation about the y 
(case B) and x axis (case C). 

The top panel of Fig. |22| shows the difference of the disc 
vertical heating, cr 2 t - <r“ 0 , for the original Aq-D2 case A 
(red), case B (blue), case C (green) and isolated (black) 
systems across the disc at the end of the simulation, t = 2 
Gyr. The heating across the disc is very similar. Only in 
case C the outer disc shows a reduced heating. The time 
evolution of the same quantity in the solar neighbourhood 
is shown in the bottom panel of Fig. [22] The difference is 
minimal during the initial 0.4 Gyr when no massive subhalo 
crosses the disc. The time evolution differs by 10km 2 s -2 
around 1 Gyr. However, the final heating rate is the same 
in all cases. Hence the impact of orientation plays a minor 
role in the vertical heating of the galactic disc. 


5 SUMMARY 

In this study we have attempted at understanding the sig¬ 
nificance of the impact of satellites on the dynamics of the 
Milky Way disc using realistic initial conditions for the DM 
substructures from the zoom-in cosmological simulations of 
Milky Way-like dark matter haloes. Such impacts are anal¬ 
ysed in the form of vertical heating of the galactic disc in the 
presence of satellite galaxies. This work benefits from four 
major assets: 1. N-body simulations which reach resolutions 
of 40 parsecs for the disc in the vertical direction, 2. ini¬ 
tial conditions for the multi-component Milky Way galaxy 
which resides in much better equilibrium and steady state 
in the isolated case compared to previous works, 3. realistic 
initial conditions for DM substructures and 4. statistical set 
of initial conditions performed for the first time. 




Figure 21. Top: The mean root-mean-square of the disc height 
for all 7 systems with satellites (solid), with satellites except Aq- 
F2 (dotted) and all 7 isolated system (dashed) at 2 Gyr. Bottom: 
The mean difference of the disc vertical heating between t = 2 and 
0 Gyr for similar systems as in the top panel. The shaded regions 
correspond to the lcr dispersion around the mean value for all 7 
systems with satellites (blue), with satellites except Aq-F2 (pink) 
and isolated systems (green). Filled and open diamonds represent 
the observed value of vertical heating in the solar neighbourhood 
and the observed value with 0.5 Gyr delay. 


To have a statistically robust understanding of the satel¬ 
lites infall, we used a distribution of subhaloes from seven 
cosmological simulations including six Aquarius simulations 
(at level 2 resolution) together with the Via Lactea II. We 
used the z = 0 snapshots containing a few hundred thou¬ 
sands DM sublialoes. The masses, velocities and positions 
of all the sublialoes were normalised using the ratio between 
the host DM halo masses (M 200 ) of the corresponding simu¬ 
lations and Aq-D2. The subhaloes at the low mass end of the 
distribution only consist of a few hundred particles, there¬ 
fore we decided to consider only objects with Mud ^ 10® Mg 
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Figure 22. Top: The difference of the disc vertical heating for 
original Aq-D2 case A (red), case B (blue), case C (green) and 
isolated (dotted) systems across the disc at the final snapshot, 
t = 2 Gyr. Bottom: the time evolution of the difference of the disc 
vertical heating for similar systems as in the top panel in the solar 
neighbourhood. 


for all our analysis. The fraction of total DM mass residing 
within substructures, f EU b, is a measure of the dumpiness 
and ranges from 5% for VLII to 14% within rso of the host 
halo. Such a range is an indication for the cosmic scatter be¬ 
tween the DM mass distribution of different simulations. In 
all simulations this fraction drops by three orders of magni¬ 
tude as we move towards the inner parts of the main haloes. 
The distribution of subhaloes within the zoom-in cosmolog¬ 
ical simulations possesses differences between various simu¬ 
lation suites. VLII has the least isotropic distribution of DM 
substructures compared to Aquarius. The subhaloes appear 
to be preferentially located around the 2 -axis in VLII while 
we also observe an anisotropic distribution in the x-y plane. 

The orbits of these subhaloes were integrated in the 
presence of their background host halo potential to identify 
the subhaloes which come closer than 25 kpc to the centre. 


The percentage fraction of these crossed subhaloes ranges 
from 6.9% to 14.5%. The cumulative subhalo mass distribu¬ 
tion and the maximum circular velocity distributions of the 
crossed subhaloes show a similar slope as the total popula¬ 
tion, in the case of all 7 simulations, implying the crossed 
sample is a good representative of the total. The crossed can¬ 
didates located at r > 100 kpc contribute to < 20% of total 
number of subhaloes at the corresponding radii, while more 
than 50% of the total crossed objects originate from this 
outer region. There exists an increasing trend for the con¬ 
centration of subhaloes, ctid = rtid/r s , with increasing dis¬ 
tance. Crossed candidates located at r > 100 kpc have lower 
c t id than the total population. Considering the large con¬ 
tribution of objects at this distance to the crossed sample 
together with higher mean concentration of the total case, 
crossed candidates have a shift towards lower Ctid across the 
whole mass range by a factor of 1.35. 


The impact strength of these subhaloes is expected to 
scale oc M 2 id , hence we chose only the crossed subhaloes with 
A/tid ^ 10 8 Mg to use in our N-body simulations. Aq-F2 has 
the most number of candidates satisfying this criteria with 
24 sublialoes, while Aq-E2 has the least with only 12. Four 
subhaloes with M t id ^ 10 9 Mq reside in Aq-D2 and Aq-F2 is 
the only simulation that has one candidate with M t i d ^ 10 10 
Mq . The initial conditions for our Milky Way galaxies were 
generated using the GallC code(Yurin & Springel 20141 
with a thin exponential disc, a Hernquist bulge and a match¬ 
ing Hernquist DM halo for the original NFW host haloes. 
The parallelised N-body particle-mesh code SUPERBOX-10 
was used for performing the simulations with three nested 
grids of different resolutions per galaxy. In the case of the 
disc we used 10 million particles in order to suppress any 
unwanted perturbations caused by low number statistics in 
the disc. The isolated MW systems were evolved for 2 Gyr 
so that the components can adapt to their environment and 
also reach an equilibrium state. The subhalo candidates with 
Mtid ^ 10 8 Mq for each simulation were inserted according 
to their initial positions and velocities extracted from the 
cosmological simulations. These systems were evolved for 2 
Gyr while the isolated cases were also simulated in parallel 
as control group. The vertical heating of the galactic disc was 
determined by measuring the root-mean-square of the disc 
thickness 2 rms and the vertical velocity dispersion squared 
cr? in the solar neighbourhood and across the disc. We ob¬ 
serve a flaring in the outer parts of the disc where the mean 
value at R = 15 kpc after 2 Gyr for increase of the thickness 
due to satellites is about 100 pc compared to the isolated 
case, including all seven simulations. However, if we exclude 
Aq-F2 the increase is only 30 pc. If we take into account the 
la statistical dispersion around the mean, the thickness of 
the disc can reach 700 pc at the outer parts at 2 Gyr. It 
is important to note that it is easier to increase the verti¬ 
cal thickness of the particle distribution in the outer regions 
of the disc than in the inner parts since these particles are 
less strongly bounded to the disc. The mean value of <r 2 t - 
alfi at the end of our simulation drops strongly in the inner 
4 kpc of the disc while becoming constant as we move out¬ 
wards. This heating is 40 km 2 s -2 for all the simulations and 
< 20 km 2 s~ 2 excluding Aq-F2. This means that the specific 
energy input is homogeneously distributed over the disc in 
the radial range 4-15 kpc. The isolated system also expe- 
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riences self-heating of a similar value as the systems with 
satellites, neglecting Aq-F2. 

According to the work done by |Holmberg et al. ( |2009| ), 
the observed vertical heating in the solar neighbourhood af¬ 
ter 2 Gyr is expected to be 144 km 2 s~ 2 . The time evolution 
of the heating in the solar neighbourhood was analysed and 
a slowly increasing trend for cr 2 t - <t 2 0 was observed in all 
our systems. For the first 0.4 Gyr a minimal heating was 
observed, since no massive satellite with M t id ^ 5 x 10 8 
Mg crosses the disc within this period. However, the heat¬ 
ing rate in the initial 1 Gyr is slightly higher than, in the 
final 1 Gyr. The observed value of 144 km 2 s -2 lies more than 
3 a above the mean vertical heating taking into account all 
simulations. According to our analysis, the choice of the ori¬ 
entation of satellites distribution for the Aq-D2 case affects 
the measured vertical heating by less than 8% after 2 Gyr 
in the solar neighbourhood, which is lower than the self¬ 
heating of the isolated system. The orientation plays only a 
minor role in heating the disc vertically. 

The measured heating from this analysis is sub¬ 
dominant and significantly smaller than in previous studies 
(e.g. Ardi et ak||2003 Kazantzidis et al.|2009 l. This mainly 
arises from the fact that our cosmological simulations do 
not include many subhaloes with M t a/Mdisc > 0.2, as was 
the case in other works, and are thought to heat the disc 
more significantly. Moreover, the DM dumpiness in Aquar¬ 
ius and VLII is much smaller than that in earlier cosmo¬ 
logical simulations, as they suffered from lower mass resolu¬ 
tion. The insufficient heating in our simulations could also 
be subject to further reduction if we allow for the presence 


of gas in the disc(Moster et al. 2010). Alternative mecha¬ 


nisms are more likely to dominate the disc heating. Heating 
by transient spiral arms or by disc growth via mass accretion 


are discussed in other works (e.g. Villalobos, Kazantzidis & 


|Helmi|[2blO| |Martig, Minchev fe Flynn||2013| ). The increase 

of disc surface density results in enhanced a z with contri¬ 
bution towards heating, ranging from less than 5% to 60% 
for increasing surface density - playing an important role in 
simulations which allow for disc growth. Coherent bending 
waves and buckling instabilities are other possible scenar¬ 
ios responsible for increasing the vertical dispersion of disc 
particles although Toomre (19831 and Sellwood et al.| ( |1998| ) 
have shown such mechanisms are inefficient in terms of re¬ 
producing the observed vertical heating. Also |Lacey| ( |1984| ) 
and |Ida, Kok uba & Makino (19931 found that the scattering 
of disc stars due to GMCs is able to redirect the peculiar ve¬ 
locities of stars out of the plane however the inputted energy 
is ineffective in increasing the velocities. 

The results of this study show that the impact from 
the subhaloes on the vertical heating of the galactic disc is 
highly dependent on the high-mass end of the cosmologi¬ 
cal subhalo distribution - M t id ^ 10 u M®. The differences 
between different simulations can be thought as cosmic vari¬ 
ance calculated as the statistical dispersion around the mean 
value from all simulations. If we exclude the Aq-F2 simula¬ 
tion as the only simulation with 1 subhalo of Mad = 5.9 x 
10 10 Mg, the impact of satellites reaches only 3 km 2 s -2 at 
solar neighbourhood after 2 Gyr which is negligible. The 
best measure of high mass satellites impact is the flaring of 
the disc although no strong heating is measured. Such in¬ 
significant heating suggests high mass subhaloes are not re¬ 
sponsible for the observed heating. We conclude that 10-15% 


contribution to the observed heating may originate from the 
impact of subhaloes. The next stage of this work is an explo¬ 
ration of the remaining processes which contribute towards 
the galactic disc heating, since we now have quantified the 
share from the infall of satellites using statistically relevant 
and high resolution analysis. 
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